# This file assesses the relationship between perceptions of local communities
# and outcomes conditional on matched pair USING DIVERSITY INDEX
library(here)
library(tidyverse)
library(lme4)
library(brms)

load(here("Design", "matches_anyDA_Diversity_new.rda"), verbose = TRUE)
load(here::here("Data", "wrkdatOwnMap_anyDA_new.rda"), verbose = TRUE)

wdat0 <- inner_join(wrkdatOwnMap_new, matches_anyDA_Diversity_new)

## Add information about which unit is ranked higher versus lower within pair
wdat0 <- wdat0 %>%
  group_by(pair) %>%
  mutate(
    perc_more = rank(subjcomm.diversity.index) - 1,
    social_capital01_rank = rank(social.capital) - 1,
    community_resp01_rank = rank(community.resp01) - 1
  ) %>%
  ungroup()

## Diversity Index of own perceptions of own map
# from survey.geo.merge.R
##
## ## This next one is the diversity of the perceptions of their own maps
## survey.geo$subjcomm.diversity.index <- (survey.geo$own.community.community.percentage.chinese/100)^2+
##   (survey.geo$own.community.community.percentage.east.indian/100)^2+
##   (survey.geo$own.community.community.percentage.black/100)^2+
##   (survey.geo$own.community.community.percentage.latin/100)^2+
##   (survey.geo$own.community.community.percentage.other.asian/100)^2+
##   (survey.geo$own.community.community.percentage.aboriginal/100)^2
##
## summary(survey.geo$subjcomm.diversity.index)

## new models with new matches

baselmerfmla2 <- as.formula("~subjcomm.diversity.index+ da_prop_vm_20pct_06+(1|pair)+(1|dauid)")

lmer1di <- lmer(update(baselmerfmla2, social.capital01 ~ .), data = wdat0)
summary(lmer1di)

lmer2di <- lmer(update(baselmerfmla2, community.resp01 ~ .), data = wdat0)
summary(lmer2di)

## since was Singular above, look at fully bayesian results as a check (they are the same)

# blmer1_community_efficacy <- brm(community.resp01 ~ subjcomm.diversity.index + da_prop_vm_20pct_06 + (1 | pair) + (1 | dauid),
#   data = wdat0,
#   chains = 4, iter = 5000, cores = 4,
#   control = list(adapt_delta = 0.99, max_treedepth = 20)
# )
# summary(blmer1_community_efficacy)
# default_prior(blmer1_community_efficacy)
#
#                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# subjcomm.diversity.index    -0.13      0.02    -0.17    -0.08 1.00    19078     8412

lmer1resdi <- c(Estimate = fixef(lmer1di)[[2]], ci = confint(lmer1di, level = .95, parm = "subjcomm.diversity.index"))
lmer2resdi <- c(Estimate = fixef(lmer2di)[[2]], ci = confint(lmer2di, level = .95, parm = "subjcomm.diversity.index"))

res <- rbind(
  "Social cohesion" = lmer1resdi, # social.capital01
  "Collective efficacy" = lmer2resdi # community.resp01
)

res

save(res, file = here("Analysis", "analysis_anyDA_Diversity_new.rda"))
